function [x mVar Se Avar f condition it J]=Newton(FUN,x,varargin) % varargout
% maximises FUN, starting at x by Newton-Raphson method
tol=1e-5; maxit=100; smalleststep=.5^20;
it=1; condition=1; improvement=1; k=length(x);
%[f g H J]=feval(FUN,x,varargin{:}); %varargout
[f g H mVar Se Avar] =feval(FUN,x,varargin{:}); %varargout
while it<=maxit && condition==1 && improvement==1;
    [s1 s2]=size(H); if s1==s2 && s2>1 d=-pinv(H)*g; else d=-g./H; end      
    step=1; improvement=0;
    while step>=smalleststep && improvement==0;
        %[ff gg HH JJ]=feval(FUN,x+step*d,varargin{:}); %varargout
        [ff gg HH mmVar SSe AAvar] =feval(FUN,x+step*d,varargin{:}); %varargout
        if (ff-f)/abs(f)>=-1e-5
            improvement=1; condition=sqrt(step*step*(d'*d))>tol & (ff-f)>tol;
            x=x+step*d; f=ff; g=gg; H=HH; mVar = mmVar; Se = SSe; Avar = AAvar;
        else
            step=step/2;
        end
        %if length(x)>length(x2); return; end
    end
    it=it+1;
end
it=it-1;